#########################################################
# Name: Trevor Osaki                                    #
# Date: Oct 11, 2023                                    #
# Purpose: Test some hypotheses (Appendix 7.4)         #
#########################################################

# Here, we should be focusing on all moderates and liberals
# Everything is arranged similarly as alpha_appendix_A7_3.R 

# Set up some parameters
# These are taken from Stata
b2 = -0.453125;
w2 = -0.409396;
b1 = -0.8991416;
w1 = -0.1087866;

est = c(b2,w2,b1,w1);

v_b2 = 0.01616309;
v_w2 = 0.01674765;
v_b1 = 0.01055636;
v_w1 = 0.00991529;

v_b2_w2 = 0;
v_b2_b1 = 0;
v_b2_w1 = 0.00627924;
v_w2_b1 = 0.00793496;
v_w2_w1 = 0;
v_b1_w1 = 0;

var = c(v_b2, v_b2_w2, v_b2_b1, v_b2_w1, v_b2_w2, v_w2, v_w2_b1, v_w2_w1, v_b2_b1, v_w2_b1, v_b1,v_b1_w1,v_b2_w1, v_w2_w1,v_b1_w1,v_w1);
  
alpha = seq(-1,2,0.001);

q = 1;
k = 4;

# Set up matrices 
test2  = data.frame(matrix(0, nrow = length(alpha),ncol = 5))
colnames(test2) <- c("alpha", "wald_WB", "p_WB","wald_BW","p_BW")

theta = matrix(est,nrow = k,ncol = 1,byrow = TRUE);
r     = matrix(0,nrow = q,ncol = 1);
V     = matrix(var,nrow = k,ncol = k,byrow = TRUE);

for(i in 1:length(alpha)){
  
  test2$alpha[i] = alpha[i];
  
  scalars1 = c(1,0,-alpha[i],-(1-alpha[i]));
  scalars2 = c(0,1,-(1-alpha[i]),-alpha[i]);
  
  R1 = matrix(scalars1,nrow = q,ncol = k);
  R2 = matrix(scalars2,nrow = q,ncol = k);
  
  R1t = t(R1);
  R2t = t(R2);
  
  B1 = solve(R1%*%V%*%R1t);
  B2 = solve(R2%*%V%*%R2t);
  
  A1 = R1%*%theta - r;
  A2 = R2%*%theta - r;
  
  A1t = t(A1);
  A2t = t(A2);
  
  wald1 = (A1t%*%B1%*%A1);
  wald2 = (A2t%*%B2%*%A2);
  
  # Back out CI for W-B switchers via calculating p-value
  # The CI should be [0.155,0.791].
  test2$wald_WB[i] = wald1;
  p_wb = 1-pchisq(wald1,df = q);
  test2$p_WB[i] = p_wb;
  
  # Back out CI for B-W switchers via calculating p-value
  # The CI should be [0.348,1.033].
  test2$wald_BW[i] = wald2;
  p_bw = 1-pchisq(wald2,df = q);
  test2$p_BW[i] = p_bw;
}


# Only keep "test2"
rm(list=setdiff(ls(), "test2"))

# 1. Find the alpha that yields the highest p-value
cat("alpha=", test2[which.max(test2$p_WB),]$alpha, "for the White-to-Black switchers.") # 0.436
cat("alpha=", test2[which.max(test2$p_BW),]$alpha, "for the Black-to-White switchers.") # 0.62

# 2. Find p-values when alpha = 0, 0.5, 1
selectedAlpha2 <- test2[test2$alpha == 0.000 | test2$alpha == 1.000 | test2$alpha == 0.500, 
                      c('alpha', 'p_WB', 'p_BW')]
selectedAlpha2$p_WB <- round(selectedAlpha2$p_WB, 3)
selectedAlpha2$p_BW <- round(selectedAlpha2$p_BW, 3)
selectedAlpha2

# 3. Back out CI for W-B switchers via calculating p-value
  # This uses the fact that the Wald is distributed chi2(q).
  # The alphas with p_wb~0.05 form the approx. bounds
  # The CI should be [0.154,0.792].
p_WB_indice_low2 <- which(diff(sign(test2$p_WB - 0.05)) == 2)
p_WB_indice_up2 <- which(diff(sign(test2$p_WB - 0.05)) == -2)
p_WB_indices2 <- c(p_WB_indice_low2, p_WB_indice_low2+1, p_WB_indice_up2, p_WB_indice_up2+1)
test2[p_WB_indices2, c('alpha', 'p_WB')]

# 4. Back out CI for B-W switchers via calculating p-value
  # The alphas with p_bw~0.05 form the approx. bounds
  # The CI should be [0.347,1.034]
p_BW_indice_low2 <- which(diff(sign(test2$p_BW - 0.05)) == 2)
p_BW_indice_up2 <- which(diff(sign(test2$p_BW - 0.05)) == -2)
p_BW_indices2 <- c(p_BW_indice_low2, p_BW_indice_low2+1, p_BW_indice_up2, p_BW_indice_up2+1)
test2[p_BW_indices2, c('alpha', 'p_BW')]
